%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,

%  the interaction matrix is set up and diagonalized as a whole.

%

% Arbitrary spin structure and comparison with Breit-Rabi formula.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function Breit

Constant;

figure(1);							clf;

B=[0:.1:5]'*1e-4;



% I =1;		L =0;		J =1/2;		S =1/2;		Aj=152.137e6;			Bj=-.1;						%6Li 2S1/2

% %E1=BreitRabiFormel(B,I,J,L,S,Aj,Bj);

% E =PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj);

% 

% subplot(221),plot(1e4*B,1e-6*E','-k','LineWidth',1.3)%,1e4*B,1e-6*E1,'-r'

% set(gca,'xLim',[0 max(B*1e4)],'FontSize',12,'LineWidth',1.3);

% xlabel('\fontsize{12}{\itB}   (G)')

% ylabel('\fontsize{12}\Delta{\it\nu}   (MHz)')





I =3/2;		L =1;		J =3/2;		S =1/2;		Aj=-(I+.5)/2*1.17e6;		Bj=0;						%6Li 2P3/2

E =PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj);



subplot(222),plot(1e4*B,1e-6*E','-k','LineWidth',1.3)

set(gca,'xLim',[0 max(B*1e4)],'FontSize',16,'LineWidth',1.3);

xlabel('\fontsize{16}{\itB}   (G)')

ylabel('\fontsize{16}\Delta{\it\nu}   (MHz)')





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,

%  this routine uses the Breit-Rabi formula. The labeling is define by

%

% mImJ= 1  2   .. 2I+1 2I+2 .. 2I+1+2J+1

% mJ  =-J -J   .. -J   -J+1 .. +J

% mI  =-I -I+1 .. +I   -I	 .. +I



function vG = BreitRabiFormel(B,I,J,L,S,Aj,Bj)



Constant;

vHFS=2*Aj;

gJ=2;

F =1;

mF=[F:-1:-F];						[mF_,B_]=meshgrid(mF,B);

x=uB*gJ*B_/(2*pi*hP)/vHFS;			vG1=-vHFS/2/(2*I+1) - vHFS/2*sqrt(1+4*mF_/(2*I+1).*x+x.^2);

vG1=vG1;



F =2;

mF=[F:-1:-F];						[mF_,B_]=meshgrid(mF,B);

x=uB*gJ*B_/(2*pi*hP)/vHFS;			vG2=-vHFS/2/(2*I+1) + vHFS/2*sqrt(1+4*mF_/(2*I+1).*x+x.^2);

vG2=vG2;

vG=[vG1 vG2];





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,

%  this routine settles the matrix and gets its eigenvalues.



function E = PaschenBackHFSMatrix(B,I,J,L,S,Aj,Bj)



for k=1:length(B)

	for mImJ1=1:(2*I+1)*(2*J+1)

		for mImJ2=1:(2*I+1)*(2*J+1)

%			mJ1=floor((mImJ1-1)/(2*I+1))-J;			mI1=mImJ1-I-1-(J+mJ1)*(2*I+1);

%			mJ2=floor((mImJ2-1)/(2*I+1))-J;			mI2=mImJ2-I-1-(J+mJ2)*(2*I+1);

			mI1=floor((mImJ1-1)/(2*J+1))-I;			mJ1=mImJ1-J-1-(I+mI1)*(2*J+1);

			mI2=floor((mImJ2-1)/(2*J+1))-I;			mJ2=mImJ2-J-1-(I+mI2)*(2*J+1);

			Mhf(mImJ1,mImJ2)=PaschenBackHFSCoefficient(B(k),I,J,L,S,mI1,mJ1,mI2,mJ2,Aj,Bj);

		end

	end

	E(:,k)=sort(eig(Mhf));

end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Zeeman and Paschen-Back splitting of a fine structure level and shift of the hyperfine structure sublevels,

%  this routine calculates one matrix element.



function Chf = PaschenBackHFSCoefficient(B,I,J,L,S,mI1,mJ1,mI2,mJ2,Aj,Bj)



Constant;

gJ=(3*J*(J+1)+S*(S+1)-L*(L+1))/2/J/(J+1);		gI=gJ*uK/uB;



denom=2*I*(2*I-1)*2*J*(2*J-1)+1e-9;

Im =sqrt(I*(I+1)- mI2   *(mI2-1));				Jm =sqrt(J*(J+1)- mJ2   *(mJ2-1));

Ip =sqrt(I*(I+1)- mI2   *(mI2+1));				Jp =sqrt(J*(J+1)- mJ2   *(mJ2+1));

Ipp=sqrt(I*(I+1)-(mI2+1)*(mI2+2));				Jpp=sqrt(J*(J+1)-(mJ2+1)*(mJ2+2));

Imm=sqrt(I*(I+1)-(mI2-1)*(mI2-2));				Jmm=sqrt(J*(J+1)-(mJ2-1)*(mJ2-2));

if mI1==mI2       & mJ1==mJ2

	Chf=(Aj+3*Bj/denom)*mI2*mJ2 - Bj*2*I*(I+1)*J*(J+1)/denom + (gJ*uB*mJ2-gJ*uK*mI1)*B/(2*pi*hP) + 6*Bj/denom*((mI2*mJ2)^2+(Im*Jp)^2/4+(Ip*Jm)^2/4);

elseif mI1==mI2+1 & mJ1==mJ2-1

	Chf=(Aj+3*Bj/denom)/2*Ip*Jm + 6*Bj/denom*Ip*Jm*(mI1*mJ1+mI2*mJ2)/2;

elseif mI1==mI2-1 & mJ1==mJ2+1

	Chf=(Aj+3*Bj/denom)/2*Im*Jp + 6*Bj/denom*Im*Jp*(mI1*mJ1+mI2*mJ2)/2;

elseif mI1==mI2+2 & mJ1==mJ2-2

	Chf=6*Bj/denom*Ip*Jm*Ipp*Jmm/4;

elseif mI1==mI2-2 & mJ1==mJ2+2

	Chf=6*Bj/denom*Im*Jp*Imm*Jpp/4;

else

	Chf=0;

end